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ABSTRACT 

Simulations of the magnetorotational instability (MRI) in 'unstratified' shearing boxes 
exhibit powerful coherent flows, whereby the fluid vertically splits into countcrmov- 
ing planar jets or 'channels'. Channel flows correspond to certain axisymmetric linear 
MRI modes, and their preponderance follows from the remarkable fact that they are 
approximate nonlinear solutions of the MHD equations in the limit of weak mag- 
netic fields. We show in this paper, analytically and with one-dimensional numerical 
simulations, that this property is also shared by certain axisymmetric MRI modes 
in vertically- stratified shearing boxes. These channel flows rapidly capture significant 
amounts of magnetic and kinetic energy, and thus are vulnerable to secondary shear 
instabilities. We examine these parasites in the vertically stratified context, and esti- 
mate the maximum amplitudes that channels attain before they are destroyed. These 
estimates suggest that a dominant channel flow will usually drive the disk's magnetic 
field to thermal strengths. The prominence of these flows and their destruction place 
enormous demands on simulations, but channels in their initial stages also offer a 
useful check on numerical codes. These benchmarks are especially valuable given the 
increasing interest in the saturation of the stratified MRI. Lastly we speculate on 
the potential connection between 'run-away' channel flows and outburst behaviour in 
protostellar and dwarf nova disks. 
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1 INTRODUCTION 

Differential rotation and weak magnetic fields conspire to 
destabilise laminar flow in accretion disks (and other as- 
trophysical settings) via the magnetorotational instability 
(MRI) (Balbus and Hawley 1991). The MHD turbulence 
that this linear instability sustains can transport significant 
angular momentum, and thus drive the mass accretion that 
is observed in many disk systems (Hawley et al. 1995, Stone 
at al. 1996, Balbus and Hawley 1998, Hawley 2000). 

The characteristics of MRI turbulence have been stud- 
ied extensively with local three-dimensional simulations in 
the 'unstratified' shearing box 1 , which is a model that mim- 
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1 Throughout the paper we adhere to the common (mis-) use of 
the term 'stratified', by which we mean a background variation 
in density and pressure but not necessarily in entropy. Normally, 



ics a small block of gas deeply embedded in the disk, igno- 
rant of its large-scale vertical and radial structure (Balbus 
and Hawley 1998). When a net vertical magnetic flux passes 
through such a box, the fluid can be dominated by coher- 
ent structures (often recurring) called channel flows. Akin to 
Kolmogorov shear flow, these consist of several countermov- 
ing planar streams of fluid situated at different vertical levels 
(Meshalkin and Sinai 1961, Hawley et al. 1995, Sano and In- 
utsuka 2001, Sano 2007, Goodman and Xu 1994, hereafter 
GX94). Previous work has investigated the role these struc- 
tures, and their destruction, play in controlling MRI satu- 
ration in the shearing box model (GX94, Latter et al. 2009, 
hereafter LLB09, Pessah and Goodman 2009, Pessah 2009). 
In this paper we generalise these analyses to the case of a 
shearing box in which the vertical structure of the disk has 
been incorporated. 



however, 'stratified' explicitly means that the background entropy 
does vary. 



The preponderance of channel behaviour in unstratificd 
boxes is connected to the fact that the channels correspond 
not only to linear MRI modes but also to approximate non- 
linear solutions of the MHD equations, when in the incom- 
pressible, weak field limit (GX94). In this paper, we show ex- 
plicitly that analogous linear modes in a vertically-stratified 
box possess the exact same property — they are also ap- 
proximate nonlinear solutions 2 . This property has impor- 
tant consequences. It means that channel flows also dom- 
inate stratified simulations with a net vertical flux, grow- 
ing exponentially and free of nonlinear interactions for a 
time sufficient to accumulate significant amounts of kinetic 
and magnetic energy. But the ultimate destruction of such 
a large amplitude channel can be so violent that the com- 
putational domain is evacuated of mass, and the simulation 
compromised (Stone et al. 1996, Miller and Stone 2000). 
In contrast, behaviour elicited by other magnetic configura- 
tions (a purely azimuthal field, for example) is not nearly 
so striking, nor numerically problematic (Miller and Stone 
2000). As a result, published simulations with a net vertical 
field are comparatively rare, even if the magnetic configura- 
tion itself is of astrophysical importance. 

The agents responsible for the disruption of isolated 
channel flows are nonaxisymmetric parasitic instabilities, 
which feed on their velocity and magnetic shear (GX94, 
LLB09, Pessah and Goodman 2009, Pessah 2009). In gen- 
eral, these parasitic modes are ineffectual until large channel 
amplitudes are attained because of the combined action of 
the background shear and the exponential growth of the 
channel itself. In this paper, we compute directly some of 
the ideal MHD parasitic modes that might be present in 
stratified boxes in a convenient intermediate limit. Using 
these results, we estimate the minimum channel amplitude 
that permits a successful parasitic attack. In a later article 
a general survey with three-dimensional simulations will be 
presented. 

The paper is organised as follows. First, in Section 2 we 
present linear MRI channel modes in a vertically-stratified 
box and show that they are also nonlinear solutions in the 
limit of weak magnetic fields. These are directly computed 
numerically and then checked with one-dimensional strati- 
fied simulations in Section 3. Following that, in Section 4, 
the ideal MHD parasitic instabilities which could disrupt 
these channels are characterised and discussed. Our conclu- 
sions are drawn in Section 5, where we also briefly discuss 
the potential role of these solutions in real accretion disks. 



MRI is local, and as such is relatively impervious to the 
boundary conditions and larger-scale properties of the disk. 
Indeed, the nature of the turbulence induced by the MRI 
has been studied profitably in the local setting of the shear- 
ing box for over a decade. However, the inclusion of vertical 
structure can help us address important questions such as 
the relationship between the MRI, on the one hand, and on 
the other: magnetic buoyancy, depth-dependent ionisation, 
realistic radiation effects, and the formation of a disk-corona 
and large-scale magnetic fields (see for example, Branden- 
burg et al. 1995, Stone et al. 1996, Johansen and Levin 2008, 
Davis et al. 2010, Shi et al. 2010, Sano and Miyama 1999, 
Stone and Fleming 2002, Miller and Stone 2000, Hirose et 
al. 2006, Blackman and Pessah 2009). 

The first linear axisymmetric stability analysis of the 
ideal MRI in a vertically-stratified local model was under- 
taken by Gammie and Balbus (1994) (see also the full global 
treatment of Papaloizou and Szuszkiewicz 1992). They ad- 
dressed different equilibrium magnetic field configurations 
in addition to various vertical boundary conditions. Gener- 
ally, the MRI eigenmodes must be calculated numerically, 
though some analytic progress can be made when the field 
is purely vertical and if it is assumed that the gas extends 
over all z (the 'infinite disk' model). These results were gen- 
eralised in a more realistic geometry by Ogilvie (1998). A 
number of papers since have numerically computed the lin- 
ear eigenmodes, usually in the context of the non-ideal MHD 
prevalent in protostellar disks (see, for example, Sano and 
Miyama 1999, Salmeron and Wardle 2003, 2005, 2008). More 
recently, Liverts and Mond (2009) returned to the isother- 
mal ideal MHD infinite disk and derived analytic approxi- 
mations to the linear modes with matched asymptotics in 
the limit of large plasma beta. 

In this paper we confine ourselves to the simple case 
of ideal MHD and an infinite disk threaded with a vertical 
magnetic field. Though the axisymmetric linear MRI modes 
in this setting have been well-established, what has perhaps 
not been fully appreciated is that for large midplane /3 a sub- 
set of these linear modes (the channel flows) are approximate 
nonlinear solutions as well, a fact independent of the partic- 
ulars of the vertical stratification. This property of stratified 
MRI channel flows has been suggested by some simulations 
(Miller and Stone 2000) and by 'hybrid' MRI-wind solu- 
tions (Ogilvie and Livio 2000, Salmeron et al. 2007). Here 
we demonstrate the result explicitly. 



2 MRI CHANNEL MODES AS NONLINEAR 
SOLUTIONS 

We study the MRI in a 'quasi-global' setting, in which the 
vertical structure of the disk is included but the radial struc- 
ture is neglected (except, of course, for the Keplerian shear). 
In the presence of weak magnetic fields, the essence of the 

2 Alternatively one could state that the Goodman and Xu chan- 
nels are exact nonlinear solutions of the equations of incompress- 
ible MHD, and the solutions we present in this paper are exact 
nonlinear solutions of the equations of anelastic MHD. But, since 
both classes of solution leave each regime so rapidly, we prefer to 
call them approximate solutions to the full equations. 



2.1 Governing equations 

We study a disk of perfectly conducting fluid orbiting a mas- 
sive central point with a frequency $l(r) (where r is cylindri- 
cal radius). Consequently, we employ the equations of ideal 
MHD and frame these in the geometry of the stratified shear- 
ing sheet (Goldreich and Lynden-Bell 1965, Gammie and 
Balbus 1994). This is the standard local approximation that 
describes a 'block of disk' centred on the midplane and at 
a radius ro, and moving along the circular orbit that radius 
prescribes. The Cartesian coordinates x and y of the block 
correspond to the radial and azimuthal directions respec- 
tively. The differential rotation is accounted for by the Cori- 
olis force and a background linear shear, u — —2Aoxe y , 



where 
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A Keplerian disk yields Ao = 3/4 and throughout this paper 
a Keplerian rotation law will always be assumed. In addition, 
to account for the vertical gravity of the central object, we 
include the force —fl 2 z e z in the momentum equation (with 
fio = Q(ro)), which is simply the first term in an expansion 
of the gravitational potential in z. The approximation breaks 
down for large z, therefore we have assumed that the disk is 
relatively thin. From now the subscript '0' will be dropped. 
The governing equations of the problem are 



dtp = -V • (pit), 

d t u + u ■ Vit = — V<E> — 2Q,e z x u 



(1) 
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d t B = Vx(uxB), (3) 
where the tidal potential is defined through 



d> = -fnV + ifiV. 



(4) 



To this set we must add the equation of state and possibly 
the internal energy equation. We must also ensure that the 
magnetic field B is solenoidal. 



2.2 Equilibrium and solution ansatz 

Let us suppose the disk is pierced by a weak vertical mag- 
netic field. At some large height the magnetic field lines will 
most likely arc back and reenter the disk at a different ra- 
dius, or perhaps plunge into the central object itself, but we 
neglect this larger structure, and consequently ignore any 
local force that a magnetic field line may exert due to the 
influence of its remote 'footpoint'. We also ignore the pos- 
sibility that the field lines are embedded in a hot magnetic 
halo. (Both cases are treated in Gammie and Balbus 1994.) 
The governing equations then admit the following equilib- 
rium state: 



uo 



(3/2)Qxe y , 
ft = ft (|). 



pa = poo h \Jj 
Bq — B e z 



In the above, h is a dimensionless symmetric function and H 
is the characteristic lengthscale of the disk's vertical struc- 
ture. The dimensional poo and Bo are constants. In addition, 
there may be variations in temperature and entropy. The two 
functions h and Po can be determined from the steady force 
and energy balances once the equation of state is specified, 
but their exact forms are unimportant in what follows. We 
require only that the stratification is convectively stable. 

Consider a perturbation to the basic state, dependent 
on only z and t, and taking the form of a planar channel 
flow: 

it ch = buo F (^jj^j e s * l e x cos# + e y sin 6], (5) 
B ch = bB oG(jj) e s * [e x sin (9- e y cos 6], (6) 
where F and G are dimensionless functions, s is a growth 



rate, 9 is a constant orientation angle, and «o and Bo con- 
stants, the latter introduced earlier in the equilibrium solu- 
tion. Note also the dimensionless constant b, which measures 
the relative amplitude of the channel to its background at 
t = 0. These perturbations constitute the basic MRI. 

It is easy to check that the profiles Eqs (5)-(6) possess 
the following attractive properties: 



ch T-T ch ?-jch 
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VB ch = 0, 



0. 



These cancellations mean that all the nonlinear terms in 
the induction equation and most of the nonlinearities in the 
momentum equation disappear. 

Via magnetic pressure, the MRI channel flow drives per- 
turbations in the thermodynamic variables, p ch , P ch , etc, 
but these will remain small relative to the background equi- 
librium for as long as the magnetic field (the driving) is 
weak. We assume that the channel's magnetic field is indeed 
weak, at least initially, and consequently linearise the gov- 
erning equations in the thermodynamic perturbations, more 
specifically p ch . The continuity equation becomes 

chx 



V- (hu 



0. 



(7) 



to leading order, which is satisfied by construction, and the 
thermodynamic perturbations drop out entirely from the x 
and y components of the momentum equation, i.e. 1/p is 
approximated by 1/po in the magnetic tension terms and 
we obtain 
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(Now all nonlinearities have vanished from these equations.) 
As a consequence, the z component of the momentum equa- 
tion, 
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the equation of state, and the entropy equation are sub- 
ordinate to the MRI dynamics. This set of 'passive' equa- 
tions are acted upon by the MRI equations, via the channels' 
magnetic pressure, but not vice-versa. Their sole job is to 
determine the small thermodynamic fluctuations driven by 
the channels and we are not obliged to solve them explicitly. 
Our only concern is for the magnitudes of these fluctuations, 
which must remain small or else the approximation breaks 
down. From (9) it is clear that these amplitudes, relative to 
their background state, are of order e — b 2 e 2st //3, where 

^Wi^oVy (10) 

is the (midplane) plasma beta associated with the back- 
ground vertical field. Our linearisation in p ch is only valid 
for as long as e remains small. 

This formalism is actually the anelastic approximation, 
modified so as to accomodate the MRI (Gough 1969, Gilman 
and Glatzmaier 1981, Barranco and Marcus 2005). Notice 
that the incompressible limit employed by GX94 is a special 
case, obtainable if we further assume the channel flows vary 
on scales much smaller than H . 

We conclude that the channel flow ansatz introduced 
in Eqs (5)- (6) is not only a linear solution but also an ap- 
proximate nonlinear solution to the governing equations in 



the regime of weak magnetic field (or, equivalently, small 
thermodynamic fluctuations). Because the approximation 
introduces errors of order e 2 , we expect the ansatz to fail 
when e ~ 1. Owing to the channels' exponential growth, 
this regime is achieved rapidly unless other agents intervene 
(such as parasitic instabilities). The time at which the ap- 
proximation fails can be simply estimated by 



t 



2s 



ln(/3/6 2 



(11) 



With a typical MRI growth rate and b = 0.1 and /3 = 10 4 , 
this time corresponds to only 1.5 orbits. But it also means 
that, in principle, a single MRI stratified channel can gen- 
erate equipartition strength fields from very weak fields in 
astonishingly short times, just as its unstratified counter- 
parts can (GX94, Hawley et al. 1995, LLB09). This is in fact 
witnessed in the stratified simulations of Miller and Stone 
(2000). 



2.3 The eigenproblem for channel flows 

We now supply the remaining details of the ansatz by de- 
riving equations for F and G and a dispersion relation for s. 
The solution (5)-(6) is substituted into the governing equa- 
tions so that u = uq + u cb and B = B + B ch . After lin- 
earising in p ch , the four equations governing the MRI chan- 
nel flow are the x and y components of the momentum and 
induction equations. These are 

u (s cos (9 - 2Q sin (9) F = (v A ) 2 s ^°\T^-j , (12) 

V ft LIZ J 

w o (ssin0 + ificos0)F = -(u A ) 2 cos<9 \\^- \ , (13) 
2 \h dz I 



■an o dF 

s sin G ■ = uo cos — — , 

dz 

3 dF 
(— s cos + — sin 5) G = «o sin# — — , 



(14) 
(15) 



where va = Bo/y/Airpoo is the equilibrium Alfven speed at 
the midplane. The four functional equations can be satisfied 
if 

f- 1 dG r- 1 dF (m) 

which together return the linear second-order ordinary dif- 
ferential equation: 



*L+K*hF = Q, 
dz 2 



(17) 



where the 'vertical wavenumber' K, a constant, is yet to be 
determined. Note that Equation (17) was also derived by 
Gammie and Balbus (1994) (cf. their Eq. (23)). 

Equation (17) with suitable boundary conditions de- 
scribes a ID eigenvalue problem for the vertical MRI 
modes with K as eigenvalue. Because h is always posi- 
tive the equation is in classical Sturm-Liouville form and 
so we are assured of a discrete set of real eigenvalues 
{K n } and eigenfunctions {F n }. Moreover, the eigenfunc- 
tions are orthogonal under integration with weight h (Arfken 
1970). The solutions possess the symmetry (K n , F n , G„) — > 
(—K n , F n , —G„). Notice that the lowest order solution is 
trivial, Ko = and F a constant, but this does not corre- 
spond to the MRI. 



The boundary condition comes from the requirement 
that the Alfven speed of the perturbation (v'a) 2 oc G 2 /h 
must never diverge. Thus G (and hence dF/dz) must go 
to zero at the disk boundaries. Solutions to Eq. (17) are 
computed in the following subsection for various equilibria. 

Returning to Eqs (12)- (15), the z functions factor out 
and we are left with the unstratified incompresible, but dis- 
crete, dispersion relation of the MRI. Eliminating and uo 
we get for a given n mode 

s 4 + (O 2 + 2v 2 A K 2 n )s 2 + v\Kl(v 2 A K 2 n - 3fi 2 ) = 0. (18) 

If time and space are scaled by Q~ x and H, respectively, the 
dispersion relation only depends on the two dimensionless 
quantities j3 and n. Evidently, the basic physics of the MRI 
is shared at each separate z, and the instability is not dra- 
matically altered by the stratification. Note that Eq. (18) 
was recovered in Gammie and Balbus (1994) and Liverts 
and Mond (2009). In the latter it appears in their Eq. (19), 
if one sets K n = n(n + l/2)/<6(— zo) in their notation. To 
complete the description we need and uo, which can be 
computed from 



= i 



3Q . 2 

iin = sin 

2K n 



1 + 32- 



(19) 



(20) 



The quantity tto / (HQ.) may be regarded as the Mach num- 
ber of the channel. Exactly as in the unstratified MRI, the 
fastest growing mode possesses the orientation angle near- 
est 7r/4, and the flow speed is always sub-Alfvenic (GX94, 
LLB09). Finally, Eq. (19) yields a critical f3 below which 
there can be no MRI at all: we find the critical value to be 
Pc = (32/15)(-Ki#) 2 . Because the longest mode (n = 1) 
will generally possess a wavenumber K\ > 1/H, we obtain 
/3 C > 2. This is the well-known linear result that instability 
is suppressed when the magnetic field (in fact, the magnetic 
torsional stress) is too strong (Balbus and Hawley 1991). 
Of course, this low /3 regime is not relevant for our large /3 
nonlinear solutions. 

2.4 Examples 

2.4-1 Unstratified disk 

The unstratified theory (Balbus and Hawley 1991, GX94) is 
recovered by setting h = 1 and sending the vertical bound- 
ary conditions to positive and negative infinity. Equation 
(17) is then the harmonic oscillator equation, and 



F = sin(Kz). 



(21) 



Because the boundary conditions have been dropped, the 
vertical wavenumber A" is a free parameter, and because H 
has disappeared from the problem the plasma beta has as 
well. Equation (21) also corresponds to the 'local' (WKBJ) 
solution of (17), which assumes that the eigenfunction varies 
on scales much shorter than H. 



2.4-2 Isothermal disk 

The isothermal ideal gas equation of state is P — c a p, where 
c s is the constant sound speed. This yields a Gaussian equi- 
librium density stratification: h — exp[— z 2 /(2H 2 )], with 



2.4-3 Approximate isothermal model 




Figure 1. Eigcnfunctions F n and G n of order n = 1, 2, 3, 4 and 10 
computed numerically for the isothermal model of Section 2.4.2. 
The functions are plotted with solid lines and are normalised so 
that the maximum of | G n | is 1 . The corresponding wavenumbers 
are K„H = 1.1584, 2.0796, 2.9829, 3.8798 and K 10 H = 9.2239. 
The dashed lines are the corresponding approximations of Sec- 
tion 2.4.3, using the Legendre polynomials. The corresponding 
eigenvalues are 1.4142, 2.4495, 3.4641, 4.4641, 10.4881. 



H — c s /Q. In this case the eigenvalue equation cannot be 
solved analytically. But it is a straightforward task to com- 
pute the eigenfunctions using a pseudo-spectral method. We 
solve for G rather than F as it decays as \z\ — > oo and so we 
can employ a basis of Whitakker cardinal functions (Boyd 
2001, also see Section 4.4) . MRI channel eigenfunctions are 
plotted in Fig. 1 for n = 1 to 4 and n = 10. 

As is clear, the MRI flow resemble a vertical sequence of 
planar 'channels' or 'jets'. The centre of each jet corresponds 
to a magnetic null surface. This is a general result which 
follows from the second equation in (16): a maximum of 
F is a, zero of G. The fluid velocities (F ) go to constant 
values as z — y oo, while the magnetic fields (G) go to zero, 
in accordance with the boundary conditions. When n is an 
even integer there exist n — 1 jets, one centred at z = 0, and 
symmetric pairs on each side at ±Zi for i — 1, ... (n — 2)/2. 
The Zi are the zeros of the function G n . When n is odd 
there are n — 1 jets occuring in antisymmetric pairs, and 
these are superimposed on a large-scale shear flow. In Fig. 1 
this pattern can be observed; the somewhat special n = 1 
case corresponds to a simple shear flow (no jets) , the n — 2 
case corresponds to a single jet centred at the midplane, the 
n = 3 case corresponds to a double jet superimposed on a 
shear flow, the n = 4 case to a triple jet, and so on. Finally, 
note that because K\ = 1.158, we have /3 C = 2.86. 



1 The numerical script that solves the isothermal eigenproblem 
is available on email request 



A rough approximation to the isothermal stratification is 
the profile h = sech 2 (z/.H") (which also describes a purely 
self-gravitating layer — Gammie and Balbus 1994). The ad- 
vantage of this choice is that it furnishes us with an analytic 
solution to Eq. (17). After some straightforward analysis one 
obtains 



Fn=Pn[tailh{z/H)], K n 



(22) 



where P n is a Legendre polynomial of order n (Arfken 
1970). These approximations are plotted in Fig. 1 as dashed 
lines alongside the numerical isothermal eigenfunctions. As 
is clear, the Legendre polynomials do an adequate job at 
replicating the isothermal channels, with the most notice- 
able discrepancies at large \z\. 



2.4.4 Polytropic disk 

A polytropic disk can be described by P — C p 1+1 / m where 
m is the polytropic index and C is a dimensional constant. 
Such a disk will settle into an equilibrium density stratifica- 
tion of the form h — (l — z 2 / H 2 ) m , where H is a combination 
of m, C, fi, and po, and denotes the actual disk surface — 
there exists only vacuum where \z\ > H . For general m, the 
boundary value problem (17) must be solved numerically, 
though an analytical solution in terms of parabolic cylin- 
der functions is available for the special case of m — 1. In 
Fig. 2 numerically computed solutions are presented when 
m = 3/2, which corresponds to an ideal gas with an adia- 
batic index of 5/3. These are produced by a simple shooting 
method. We plot the first four nontrivial eigenfunctions as 
well as the n = 10 solution. The critical plasma beta, be- 
low which there exists no unstable MRI, is /3 C = 13.8 when 
m = 3/2. 



2.5 Discussion 

The 'nonlinear property' of the channel modes permits the 
accumulation of significant kinetic and magnetic energies 
in the first few orbits of the MRI evolution. At large mid- 
plane P, small initial conditions seed a suite of fast-growing 
channel modes, which continue growing independently even 
once they reach the magnitude of the background vertical 
field Bo • Because channels modes do not nonlinearly interact 
amongst themselves there may be no check on their growth 
until equipartition amplitudes. 

There are perhaps two ways to avoid this outcome. 
First, parasitic modes, feeding off the MRI channels' strong 
velocity shear may destroy the channel prematurely (GX94, 
LLB09, Pessah and Goodman 2009, Pessah 2009). We ex- 
plore this in more detail in Section 4, but both unstratified 
and stratified simulations indicate that parasitic modes gen- 
erally cannot halt the channels' march to equipartition in the 
early stages of a run (Hawley et al. 1995, Miller and Stone 
2000, LLB09). Second, channel flows may interact nonlin- 
early with slower-growing radially dependent MRI modes or, 
in particular, transiently growing non-axisymmetric modes 
(Balbus and Hawley 1992, Hawley et al. 1995). In principle, 
if the amplitude of the initial conditions is set equal to or 
above that of the background B z , a non-axisymmetric MRI 
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Figure 3. Spacetime diagram of the gas density in a simulation 
performed with RAMSES for j3 = 100 in which a pure n = 7 
eigenmode is seeded at time t = 0. Density spikes develop between 
regions of strong magnetic field, which then lift up the thin veins 
of matter through magnetic buoyancy. 
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Figure 2. Eigenfunctions F n and G n of order n = 1,2,3 
and 10 computed numerically for the polytropic model of Sec- 
tion 2.4.4 with m = 3/2. The corresponding wavenumbers are 
K n H = 2.541, 4.762, 6.963, and 22.29 respectively. The profiles 
resemble those in Fig. 1 but are confined between z = ±H. 

mode (shearing wave) could destructively interfere with a 
channel flow, before it dominates everything else. Indeed, 
seeding large-amplitude noise is a common technique in un- 
stratified mean vertical flux simulations, as it averts the ini- 
tial channel 'spike' (Hawley et al. 1995, LLB09). This tech- 
nique should also serve in stratified simulations at interme- 
diate and lower /3. In very large (5 runs (as in Suzuki and 
Inutsuka 2009), the seeded noise is usually nonlinear from 
the very start. 

If a channel reaches equipartition, magnetic pressure 
will be sufficient to alter the channel structure, squeezing 
mass onto the horizontal sheets where the magnetic field 
is zero. Some of this behaviour is explored in the follow- 
ing section with one-dimensional numerical simulations. The 
stability of these 'extreme channel' were studied in the un- 
stratified context in LLB09, which showed they were undone 
by compressible and resistive parasitic modes. 



3 ISOTHERMAL ID SIMULATIONS 

In this section we present a set of ID numerical simulations 
using two different finite-volume 2nd order Godunov MHD 
codes. The purpose of these simulations is twofold. First, 
they will serve as a validation of the ansatz in Section 2, 
while conversely providing a useful numerical check on these 
codes. Their second aim is to determine how far into the 
nonlinear regime the linear channels can be followed before 
magnetic pressure effects become important. Obviously, an 
important limitation of the simulations is that they are ID. 
As such, they cannot capture the parasitic modes considered 



in Section 4, which should be studied with 3D numerical 
simulations. 

3.1 Numerical setup 

We used the MHD codes RAMSES (Teyssier 2002, Fromang 
et al. 2006) and NIRVANA (Ziegler 2004) to solve Eqs (1)- 
(3) along the z direction only. An isothermal equation of 
state is adopted and, accordingly, the density profile at the 
start of the simulation is set to a Gaussian. A pure verti- 
cal magnetic field is added to this hydrostatic equilibrium, 
whose strength is governed by the midplane /3. The com- 
putational domain extends in all cases from z — —5H to 
2 = +5H, which we found was the minimum size needed 
for accurate results. The vertical boundary conditions are 
as follows: the density is extrapolated in the ghost zones in 
order to satisfy the vertical hydrostatic equilibrium, outflow 
boundary conditions are applied on the velocities, and zero- 
gradient boundary conditions are applied to the tangential 
component of the magnetic field (in ID, the vertical field is 
constant in space and time). At t = 0, we either seed the 
isothermal eigenmodes calculated in Section 2.4.2 or their 
approximation in terms of Legendre polynomials, outlined 
in Section 2.4.3. In both cases, their amplitudes are such 
that the maximum value of the horizontal magnetic field is 
a fraction b of the vertical field Bo- 

3.2 The case /3 = 10 2 

We first consider the case /3 = 10 2 . Admittedly, this is a 
choice that does not best demonstrate the nonlinear prop- 
erty of the channel flows: the channels can grow only a few 
times Bo before hitting equipartition. However, lower /3 runs 
facilitate a clean comparison of the growth rates and eigen- 
functions, because the various modes are fewer and relatively 
well-spaced. 



Mode 




& num / 


^iium / 


number 




RAMSES 


NIRVANA 


p = 100 


n= 


=1 


0.2664 


0.2662 ± 0.0001 


0.2674 ± 0.0009 


rx= 


=2 


0.4307 


0.4302 ± 0.0001 


0.4276 ± 0.0001 


u- 


=3 


0.5502 


0.5495 ± 0.0002 


0.5471 ± 0.0002 


u- 


= 1 


0.6363 


0.6353 ± 0.0003 


0.6330 ± 0.0002 


u- 


=5 


0.6957 


0.6944 ± 0.0003 


0.6920 ± 0.0002 


u- 


=6 


0.7326 


0.7307 ± 0.0003 


0.7290 ± 0.0002 


u- 


=7 


0.7489 


0.7463 ± 0.0003 


0.7458 ± 0.0003 


u- 


=8 


0.7455 


0.7414 ±0.0004 


0.7427 ± 0.0004 


u- 


=9 


0.7214 


0.7153 ± 0.0004 


0.7148 ± 0.0005 


11= 


=10 


0.6744 


0.6653 ± 0.0003 


0.6618 ± 0.0003 


p = 1000 


n= 


=20 


0.7343 


0.7308 ± 0.0008 


0.7174 ± 0.0005 


u= 


=30 


0.7133 


0.6921 ± 0.0038 


0.6989 ± 0.0012 


P = 100-Legendre Approximation 


n= 


=4 


0.6363 


0.6357 ±0.0041 


0.6374 ± 0.0041 




=7 


0.7489 


0.7413 ± 0.0033 


0.7435 ± 0.0012 


n= 


= 10 


0.6744 


0.6677 ± 0.0089 


0.6737 ± 0.0064 



Table 1. Growth rates of the normal modes of order n (first 
column) obtained with the codes RAMSES (third column) and 
NIRVANA (fourth column). The second colum display the the- 
oretical growth rates for the purpose of comparison. The data 
under the heading 'Legendre Approximation' were obtained with 
runs in which the initial condition was an approximate isothermal 
cigcnfunction (cf. Section 2.4.3); all other runs were seeded with 
the correct isothermal eigenfunction (computed in Section 2.4.2). 



We set b — 0.01 and seed as initial conditions the 
isothermal eigenfunctions of Section 2.4.2 with n values be- 
tween 1 and 10. For this range of parameters, we find a 
resolution of N z == 256 gridpoints to be mandatory in order 
to properly reproduce the mode morphology. This amounts 
to about 25 cells per scaleheight. For each n, we measure 
the growth rates with both RAMSES and NIRVANA. They 
are reported in Table 1, in the third and fourth columns re- 
spectively, while the second column gives the expected the- 
oretical growth rate computed using Eqs (17) and (18) with 
the isothermal model. For all n, the agreement between the 
numerical results and the analytical prediction is excellent, 
cross-validating both results. It should be noted, though, 
that it is difficult to follow the n — 1 mode for longer than 
1.5 orbits (at which point the eigenmode has grown by only 
about an order of magnitude). Because a num is so small 
for this mode, any numerical error tends to seed higher or- 
der rapidly-growing eigenmodes that quickly overcome the 
seeded n — 1 mode. 

Finally, we comment briefly on the nonlinear evolution 
of the eigenmodes. A spacetime diagram of the density field 
p for the seeding of a n == 7 mode is presented in Fig. 3. No 
features are apparent when t < 1.3 orbits, corresponding 
to the weak-field, anelastic phase of the evolution. At near 
t = 1.3 orbits the amplitude of the channel magnetic field 
is already a few times that of the vertical field -Bo, yielding 
a midplane f3 of order 10. At later times, or rather smaller 
P, magnetic pressure fluctuations compress the density into 
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Figure 5. As in Fig. 3 but in the case n = 20 and P = 1000. 



thin layers that are then lifted up by magnetic buoyancy. 
This is similar to the results of past 3D numerical simula- 
tions of the MRI in the shearing box in the presence of a 
vertical flux (Miller and Stone 2000). 

3.3 The cases 8 == 10 3 and 8 = 10 4 

We also perform ID simulations with 8 == 10 3 and 8 — 
10 4 , both with b — 0.1. These larger 8 values offer more 
opportunity for the channel flows to confirm that they are 
indeed good approximate nonlinear solutions. 

When f3 — 10 3 , the fastest growing eigenmodes are 
those with n values between n == 20 and n == 30. To ad- 
equately describe these modes, a resolution of N z — 512 is 
required. This corresponds to more than 50 cells per scale- 
height. A lower resolution is unable to represent the high 
spatial frequencies of the profiles. Despite these large spa- 
tial resolution, we are unable to follow the growth of low 
order modes. This is again because their growth rates are 
so comparatively small: numerical errors seed higher order 
modes that grow faster and quickly distort their low order 
precursors. 

We present the results of two simulations, the first in- 
tialised with a n == 20 mode and the second with a n — 30 
mode. The growth rates of both are reported in Table 1. 
Similarly to the /3 == 100 runs, there is excellent agreement 
between the analytical and the numerical a; moreover the 
linear modes hold their spatial profiles up to large ampli- 
tudes. The latter point is illustrated in the left panel of 
Fig. 4, in which the vertical profile of B x /Bq is shown (solid 
line) at time i ===== 1 .1 orbit. The dashed line shows the same 
function at time t == 0, appropriately rescaled in order to 
match the amplitude of the eigenmode at t = 1.1. The two 
profiles are in good agreement, even though the mode has 
grown by about two orders of magnitude, and at t = 1.1 is 
nearly one order of magnitude greater than the background 
field -Bo. After t == 1.1, the midplane f3 dips below 10, and 
the subsequent evolution mirrors that of the 8 — 100 runs. 
This nonlinear behavior is illustrated in Fig. 5, a spacetime 
diagram for the density in the case n = 20. 




Figure 4. Left panel: Vertical profile of B x /Bq (solid line) at time t = 1.1 orbits for a simulation performed with RAMSES in which an 
exact normal mode of order n = 20 is seeded initially with = 10 3 . The dashed line shows the vertical profile of B x at t = 0, rescaled 
in order to match the amplitude of the field at t = 1.1. Right panel: as in the first panel, but with /3 = 10 4 , a seeded n = 68 mode, and 
after t = 1.3. The domain has been truncated to better display the fine-scale structure. 



When /3 — 10 4 , the fastest growing eigenmode is n = 68. 
Due to its fine-scale structure, a resolution of N z — 2048 
is required to describe it. This corresponds to more than 
200 cells per scaleheight (which in 3D would amount to a 
formidable computational challenge). As in the previous ex- 
amples, the numerical growth rate compares well with the 
analytic prediction: the former is 0.7398 ± 0.00018 (with 
RAMSES), while the latter is 0.7424. In the right panel of 
Fig. 4 the channel's vertical profile of B x /Bo is plotted when 
t — (dashed line) and when t = 1.3 (solid line), with the 
t — profile rescaled. Because the rapid spatial oscillations 
are difficult to see, we truncate the figures at ±2H. Note that 
by t — 1.3 the channel has grown by some 25 times the back- 
ground vertical field, but the agreement between the two 
profiles remains excellent. Discrepancies emerge first at in- 
termediate 2. This is because the (depth-dependent) plasma 
beta associated with a channel is somewhat smaller at a few 
scale heights and compressible physics distort the original 
profile in these layers before the others. Thus equipartition 
will first emerge not at the midplane but at a few scale 
heights. But this effect is only noticeable at larger n. 

3.4 Legendre polynomial approximation 

We also test the soundness of the approximation introduced 
in Section 2.4.3. To do that, we perform simulations in which 
are seeded the approximate eigenmodes of Eq. (22) , though 
still with the standard Gaussian profile for the equilibrium 
density. For /3 = 10 2 and n between 4 and 10 we find the ap- 
proximate modes hold their vertical profiles relatively well 
throughout the time evolution. This is illustrated in Fig. 6, 
which compares the vertical profile of B x /Bo obtained with 
RAMSES at time t — 1.2 orbits in simulations in which the 
exact eigenmodes are seeded (solid line) with the results ob- 
tained when the Legendre approximation is seeded (dashed 
line). The three panels, from left to right, respectively cor- 
respond to the case n = 4, n = 7, and n = 10. The n = 7 
mode is the fastest growing and it offers excellent agreement 
between both solutions. The agreement is acceptable in the 



other two cases, though not so striking. The discrepancy 
here follows simply from the smaller growth rates of these 
modes. The Legendre profiles, being approximations, possess 
small components in the 'eigendirect ions' of faster growing 
exact modes. Over time, these small perturbations, enjoying 
larger growth rates, will distort the evolution of the seeded 
approximation, resulting in the imperfect first and third pro- 
files of Fig. 6. Nevertheless, these deviations are controlled 
and fail to modify the growth rate significantly, as shown in 
the third part of Table 1 for the modes n = 4, 7 and 10. For 
values of n smaller than four and larger than ten, we found 
that such problems become important and the Legendre ap- 
proximation breaks down. This is a problem for larger j3, 
where there are more growing modes at larger n. At any 
given time, however, the set of Legendre polynomials offers 
a useful approximate spectral decomposition of any strati- 
fied flow into its component MRI modes. This point will be 
explored in more detail in a forthcoming numerical paper. 



4 PARASITIC INSTABILITIES 

Stratified channel flows will be subject to secondary 'par- 
asitic' instabilities (GX94, LLB09, Pessah and Goodman 
2009, Pessah 2009), which can feed on their strong velocity 
shear, and if resistivity is present, on their magnetic shear 
(both evident in Figs 1 and 2) . This section concerns the be- 
haviour of these secondary instabilities, characterising them 
mathematically in a convenient intermediate limit of chan- 
nel amplitude. Their role in the destruction of MRI channels 
is then discussed, as well as their importance in sustained 
MRI-induced turbulence. 

4.1 Intermediate asymptotic limit 

The parasitic analysis is complicated by the fact that the 
equilibrium state that we are perturbing varies with x and 
t (in addition to z). In shearing coordinates the x depen- 
dence can be removed, but the time dependence remains: 
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Figure 6. In each panel the solid line represents the vertical profile of the radial component of the magnetic field (solid line) after 
t = 1.2 orbits for seeded isothermal modes of n = 4 (left panel), 7 (middle panel) and 10 (right panel), all with /3 = 10 2 . The dashed lines 
represent results obtained when the Legendre approximations to these modes are seeded. 



the channel flow grows exponentially like e st , and the back- 
ground linear shear flow, issuing from Keplerian rotation, 
will 'shear out' non-axisymmetric disturbances over time. 

To get around these difficulties, we assume that the 
parasitic modes grow on a timescale much faster than the 
channel growth time and the Keplerian shear time (GX94). 
This, in fact, corresponds to an assumption concerning the 
amplitude of the channels themselves, and is equivalent to 
b^> 1. So the channels must take large amplitudes — much 
larger than the amplitude of the vertical equilbrium field 
Bq. But if we want to use the convenient channel profiles 
computed in Section 2, these amplitudes cannot be too big, 
for if they were, then the ansatz, Eqs (5)-(9), would fail to 
hold. More precisely, the midplane Alfven speed associated 
with the channel Vj± must be much bigger than the mid- 
plane Alfven speed of the vertical background va, and con- 
currently much smaller than the sound speed c. This leads 
to the restriction 

1 -C b < (23) 

Obviously, in order for this intermediate range of amplitudes 
to exist, the initial j3 must be sufficiently large. Once this as- 
sumption is made, we can use the vertical channel structures 
given by Eq. (17), while neglecting their exponential growth, 
the Keplerian shear, Coriolis force, and vertical background 
field. These simplifications transform the channel computa- 
tion into a ID eigenvalue problem in z. 

4.2 Governing equations 

We perturb the channel solution with a small incompressible 
disturbance, so that 

u = u + u ch + u, B = B + B ch + B', * = vf + 

where ^ denotes total pressure and the prime designates the 
parasites. For simplicity we have neglected the small varia- 
tions in the total pressure introduced by the channel itself. 
We have also ignored buoyancy effects associated with the 
background vertical pressure structure, which will have a 
small stabilising influence on the channel's strong velocity 
shear (Drazin and Reid 1982). Together with the incom- 
pressibility assumption, these approximations are particu- 
larly valid when K n is large, which corresponds to the regime 
of interest. The disturbances are substituted into the gov- 



erning equations and we linearise in the small parasitic dis- 
turbances. We next assume that b is large and that the time 
dependence of the parasites scales like dt ~ bQ. This allows 
the dropping of the background Keplerian shear uo , Coriolis 
force, and the exponential growth of the channel to leading 
order. Throughout most of the domain, this also means that 
the background vertical field Bo is also subdominant, except 
far from the midplane where B ch goes to zero. However, it 
can be shown that the contribution of the background verti- 
cal field in these regions makes little difference to the overall 
solution. Consequently, as far as the parasitic mode is con- 
cerned, the equilibrium is to leading order stationary, and 
we are permitted to take normal modes: 

u, B', *' oc j k °> x+ik * y+at } 

where k x and k y are (real) wavenumbers and a is a growth 
rate. Next, we adopt space and time units so that H = 1, 
and bQ — 1, while scaling the velocity, magnetic field, and 
pressure by buo, bBo, and 6 2 UqPoo respectively. This scaling 
removes the parameter b from the problem directly. Also, we 
solve for linear momentum density p' , defined through 

p' = h(z)u', (24) 

in place of the velocity. The linearised equations for p' , B' , 
and are now: 

op = —M [p' z d z u + i(fc •u)p' + (ifc + e z a z )*'] 

+ JM [#0.B + i(fc.B)B'] . ( 25 ) 

oB' = M \B' z d z u-p' z d z B/h+\p'{k-B)/h-\B' (fc • «)] , 

(26) 

= ifc ■ p + d z p' z -d z (\nh)p' z , (27) 

where k = k x e x + k y e y , the channel background is repre- 
sented by the 'hatted' vectors: 

u = F n (z) (e x cos 8 + e y sin 6) , 

B = G„(z) (e x sin (9 — e y cos 9) , 

and where we have introduced the Mach number 

M = u /(HQ) (28) 

which can be computed from Eq. (20). Equations (25)-(27) 
describe a one-dimensional eigenvalue problem with eigen- 
value a. The boundary conditions require that all variables 
go to zero at the disk's vertical boundaries. 



The inputs for the eigenproblem include: the function 
h (representing the equilibrium density stratification), the 
parameter (representing the strength of the vertical mag- 
netic field), the integer n (which fully defines the channel 
flow), and k x and k y (which define the horizontal variation 
of the parasitic mode), fn place of the last two it is conve- 
nient to use k = yjk% + k y and 6 k = axctaxi(k y /k x ). These 
are just the horizontal wavevector magnitude and its orien- 
tation with respect to the x and y axis. It is also helpful 
in some circumstances to use A = 9 — 8k, which measures 
the difference in orientation between a parasite's wavevector 
and the direction of the channel flow. All the other quanti- 
ties and functions that appear in the equations — K n , F„, 
G n , 6, M — depend solely on h, /3, and n, through Eqs (f 7), 
(19), and (20). 



4.3 Asymptotic analysis 

The full set of equations are solved numerically for general 
parameters, but some analytic expressions can be obtained 
in the limit of small k. These provide a check on the nu- 
merics and also offer some physical insight. The details of 
the analysis can be found in Appendix A; in this subsection 
we briefly summarise some of the results pertaining to the 
fastest growing parasitic mode. 

As in GX94 and LLB09, the parasites can be classed 
into two main types. There is the Type 1, Kelvin- Helmholtz 
or kink, mode which attacks every jet in the sequence of 
MRI channels, vertically buckling each in the familiar sinu- 
soidal pattern. Then there are the Type 2, or kink-pinch, 
modes, which kink some of the jets, but attack others via 
an alternate vertical squeezing and stretching, in a pattern 
something akin to a sausage (see LLB09). The latter 'hy- 
brid' modes always grow slower than the kink mode, usually 
by an order of magnitude. Because such modes will rarely 
destabilise a channel, and their asymptotic analysis is quite 
involved, we do not discuss them in this subsection. 

In the limit of very long horizontal wavenumber < 
k <C 1, the leading order features of the kink mode are rel- 
atively easy to ascertain. Its growth rate is given by 

a 2 = A n k 2 M 2 [cos 2 A — M^ 2 sin 2 A] , (29) 

where A n is a positive number dependent on the channel 
profile (see Appendix A) , and the Alfvenic Mach number is 
M A = u /v A = M^J (3/2. The latter depends on 0, but it 
will always be less than 1. The growth rate scales like k, and 
the modes will stabilise when 

\0 - 9 k \ > arctanMA. (30) 

The stability criterion above is exactly that of parasitic kink 
modes in the unstratified case (GX94, LLB09). The kink 
mode favours orientations of its wavevector near the direc- 
tion of the channel flow; the more k points away from u c , 
the less shear energy is available and the greater the stabil- 
ising magnetic tension. When the angle difference \6 — d k \ is 
larger than some number less than 7r/4 the mode stabilises. 
Conversely, the fastest growing modes are oriented paral- 
lel to the channel, 6 = 8 k , and being perpendicular to the 
magnetic field, are essentially hydrodynamic. 



4.4 Numerical solutions 

The Eqs (25)-(27) are solved numerically via a pseudospec- 
tral method similar to that employed in LLB09 but using 
Whittaker cardinal functions (as in Section 2.4.2) rather 
than Fourier modes. The Whittaker basis is convenient be- 
cause it ensures all solutions obey the decaying boundary 
conditions far from the midplane (Boyd 2001). If N functions 
and a (uniform) grid spacing of Sz are employed, the gov- 
erning equations reduce to a generalised algebraic eigenvalue 
problem of order 7(N/Sz), the spectrum (or some portion of 
the spectrum) of which may be obtained by the QZ algo- 
rithm, or an Arnoldi method (Golub and van Loan 1996). 
Setting N = 16 and Sz = 0.04 sufficed in most cases, though 
these must be improved upon if one is to resolve solutions 
when n is large (because of the channels' and their parasites' 
finer spatial structure). In order to simplify the problem we 
assume that the disk and the MRI modes can be represented 
adequately by the approximate isothermal model of Section 
2.4.3. 

As in LLB09, the eigenproblem is subject to many pa- 
rameters and exhibits a rich variety of mode behaviour. 
In particular, when n increases the number and complex- 
ity of the set of unstable modes is striking. However, we 
may always decompose this suite of instabilities into (a) a 
single 'pure' kink mode, and (b) a subset of slower grow- 
ing (and more physically complicated) kink-pinch modes. 
To ease their presentation, we restrict ourselves to hydro- 
dynamical modes only; that is, we set 8 = 6 k . In practice 
these turn out to be the fastest growing, and therefore the 
most dynamically significant. For the same reason, we con- 
centrate our attention primarily on the kink mode, and do 
not catalogue the physical intricacies of the many kink-pinch 
modes that emerge: they just grow too slowly. Finally, we 
specialise to the case /3 = 1000. This value permits only a 
marginal satisfaction of Eq. (23), but facilitates the calcula- 
tion. In any case, the general structure of the solutions arc 
relatively impervious to /3. This leaves the two parameters 
n and k. When (3 — 1000 there exist growing channel flows 
for all n < 38, with the fastest growing channel possessing 
n = 21. We present detailed results for only two cases n — 2 
and n = 10, as these summarise neatly the overall behaviour. 
In both cases, k is allowed to vary smoothly from to near 
K, at which point the parasitic modes stabilise. The results 
for the n = 21 mode are complicated, but are similar overall 
to the n = 10 case. 



4-4-1 The n — 2 MRI channel flow: a single jet 

The somewhat special n — 1 channel flow is a simple shear 
layer (see Fig. 1), and, as expected, yields only one para- 
sitic mode. This is the classical Kelvin-Helmholtz instabil- 
ity, which achieves its fastest growth when k « K/2 and 
is extinguished when k > K (for n = 1, K =1.1584). The 
n — 2 channel flow consists of a single jet centred on the mid- 
plane, and is in fact a close relative of the sech 2 « 'Bickley 
jet' (Bickley 1937). Consequently, it is subject to two par- 
asitic instabilities, the classical kink and pinch (see Drazin 
and Reid 1981, Biskamp et al. 1998, LLB09). The growth 
rates of each as a function of k are plotted in Fig. 7. As 
is typical, the growth rate of the kink mode is significantly 
larger than the pinch mode, and extends over a wider range 
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Figure 7. Growth rates of the unstable kink and pinch mode 
attacking the n = 2 MRI channel, which is a single magnetised 
jet centred at z = 0. The plasma beta /3 is set to 1000. The 
growth rates are plotted as functions of kH, while 8k = 8. The 
latter restriction means the modes are essentially hydrodynamic. 
The asymptotic estimate of a for the kink mode, cf. Eq. (29), is 
plotted with a dashed line. 
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Figure 8. The u' x , u' z and \E>' eigenfunctions of a kink mode 
attacking the n = 2 channel. Here k = 0.01K(= 0.0245), 8 k = 
8, and /3 = 1000. The solid lines indicate real part and dashed 
lines indicate imaginary part. For comparison, the black points 
represent the asymptotic eigefunction, cf. (A7). The growth rate 
of the mode is a/(bn) = 1.03 X 10" 4 . 



Figure 9. The u' x , u' z and ^' eigenfunctions of a pinch mode 
upon the n = 2 channel. Here k = Q.05K (= 0.12247), 8 k = 8, 
and /3 = 1000. The growth rate is <r/(M7) = 2.49 X 10" 5 +i 5.575 X 
10 -4 , the imaginary part of which is in reasonable agreement 
with the leading order asymptotic theory, which gives a = i5.75 X 
10 -4 . The black points represent the leading order Vl 1 ' asymptotic 
profile. 

of k, up to roughly K itself (where, for n = 2, K = 2.0796). 
For small k, the kink mode growth rate a scales as k, while 
the pinch appears to go as a ~ fc 3 . We also plot in Fig. 7 
the asymptotic prediction (29) with a dashed line. At low k 
it offers excellent agreement with the computed numerical 
values. 

Typical profiles of the two modes are presented in Figs 
8 and 9 at small k, alongside the asymptotic pressure pro- 
files. These low k modes should rarely be observed because 
they grow too slowly, but they allow a clearer demonstra- 
tion of the main physics and a clean comparison with the 
asymptotic theory. The kink mode in the first figure is easy 
to distinguish because it possesses a vertical velocity max- 
imum at the centre of the jet, z — 0. Correspondingly, the 
pressure pertubation at this point exhibits a strong vertical 
gradient. As a consequence, the mode elicits an alternating 
horizontal sequence of upward and downward motions along 
fc, which will buckle the jet in the familiar Kelvin-Helmholtz 
pattern. In contrast, the pinch mode in Fig. 9 exhibits a 
pressure maximum at the jet centre, while u' z changes sign. 
This corresponds to a horizontal sequence of squeezes and 
rarefactions as fluid is either drawn towards or repelled by 
the jet. 

Other 9k were tried, but only a very narrow range of 
orientations (around 9) permitted growing modes. As pre- 
dicted by the criterion (30), the kink mode stabilises when 
\9k — 9\ exceeds about 12 degrees. The interval is particu- 
larly small because the n = 2 channel is very sub-Alfvenic for 
/3 — 1000. Faster growing (larger n) channels yield greater 
Ma and, consequently, a greater range of unstable 9k- In 
common with its kink-pinch counterparts in an unstratified 
disk, the pinch mode grows for a broader range of 9k (GX94, 
LLB09, Pessah and Goodman 2009, Pessah 2009). 

Finally, note that these instabilities will fail to appear 
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Figure 10. Selected components of the kink mode eigenfunction, 
for parameters: n = 10, /S = 1000, and fc = 0.05 if (= 0.524), 
e k = 6. The growth rate of the mode is a/{bQ.) = 3.01 X 10~ 3 . The 
short vertical arrows in the central panel indicate the locations of 
the n = 10 MRI channels. 

in numerical simulations with short radial domains, such as 
one scale height H. Because Ki ~ 2.45/1/, the radial box 
length should be at least 2n/K2 ~ 2.57 H to capture the 
shortest horizontal modes. As n increases for given /3 this 
constraint becomes less of a problem: for n = 10 a radial 
length of H will comfortably fit a number of growing kink 
modes. 

44.2 The n = 10 MRI channel flow: 9 jets 

Larger n channels support a greater number of unstable par- 
asitic modes; this is simply because larger n channels' ex- 
hibit more complicated vertical structures and thus permit a 
wider array of destabilising configurations. The presentation 
of these many modes at various n is too lengthy and labo- 
rious a task for this paper. Instead we show results from a 
representative higher order channel, associated with n — 10. 

From Fig. 1, the n = 10 channel flow consists of a jet 
situated at the midplane z = and a symmetric set of 4 jets 
on each side. We find that these can harbour up to 9 parasitic 
modes. These include the kink mode, which grows the fastest 
and kinks every jet, and 8 kink-pinch modes, which differ 
in the number and location of their component kinks and 
pinches. Not every combination of kinks and pinches are 
supported by the system (if this was so there could be ~ 
2 9 such modes!). One important constraint is that no two 
pinches can attack adjacent jets, a limitation also evident 
in the unstratified problem (which admits no 'pinch-pinch' 
modes, LLB09). Some modes also leave a number of jets 
unmolested, with u' z both a maximum and zero at these 
locations. When k takes larger values this becomes more 
noticeable, and the fastest growing modes attack those jets 
nearest the midplane preferentially. 

We plot the eigenfunction of the 'pure' kink mode in 
Fig. 10 for small k. Again, we choose small k modes as 
they demonstrate the physics the clearest. We deal with 



Figure 11. Selected eigenfunction components of a kink-pinch 
mode, with parameters as in Fig. 10. The growth rate is a/(bQ) = 
5.84 X 10~ 4 — i3.18 X 10~ 3 . As earlier, the black vertical arows 
indicate the locations of the MRI channels. Here the central three 
channels are pinched as a whole and the remaining channels are 
kinked. 



the fastest growing modes in the following subsection. In 
the central panel we also indicate the centres of the nine 
channels by short vertical arrows. These show that the lo- 
cal extrema of the parasite's vertical velocity are located ex- 
actly at the channel centres. Thus the mode kinks or buckles 
every jet concurrently. Simultaneously, the channel centres 
correspond to those locations where the pressure perturba- 
tion possesses a large gradient, though this is less easy to 
see. Far clearer, is the large-scale antisymmetry in the pres- 
sure perturbation, especially in its real part. This means 
that in addition to the localised kinks upon each channel 
there exists a 'global' kink upon the entire set of jets as 
well. Variations in the eigenfunctions' large-scale 'envelope' 
are a feature of the more complicated high n parasites. In 
the unstratified analysis, the larger scale structure mani- 
fests simply through the Floquet factor e lk ' z (see GX94 and 
LLB09) . The maximum growth rate for the kink mode (when 
n = 10) is a = 0.01726 achieved when k = 0.56K. The mode 
stabilises when k > 0.97K. (Note that for n = 10, we have 
K = 9.2239.) 

We present one example of a kink-pinch eigenfunction 
in Fig. 11. As earlier, the small vertical arrows indicate the 
locations of the channel centres. The vertical velocity pos- 
sesses local extrema at all the centres except the jet at the 
midplane. At the midplane the vertical velocity changes sign, 
which means this mode pinches the central jet. The other 8 
jets, in contrast, are kinked (to varying degrees). The pres- 
sure pertubation tells us, on the other hand, that the mode 
exhibits a large-scale scale envelope through which the entire 
set of jets is pinched. This is also clear from the large-scale 
antisymmetry in the u' z profile. 



4-4-3 Fastest growing modes 



The maximum growth rates of the kink mode computed with 
our stratified model are significantly less than growth rates 
in the unstratified model (GX94). The unstratified parasitic 
kink mode reaches a maximum a ~ 0.2 bQ when attacking 
the fastest growing MR1 channel. In contrast, 

• when ft = I0 2 , the n = 7 MRI channel is dominant 
and its fastest growing parasite possesses a = 0.0596 bQ, at 
k x 0.525 K 7 . 

• When /3 = 10 3 , the n = 21 channel is dominant and the 
maximum a is 0.0329 at k w 0.575 /Gi- 

• And when (3 — 10 4 , the n = 68 channel dominates and 
(Jmax = 0.0192 bQ at k » 0.585 K 6 g. 

This discrepancy is something of a surprise. One may have 
expected that at sufficiently large f3 the unstratified and 
stratified results would converge. Instead there appears a 
scaling, with possibly the logarithm of f3, which diverges 
from the unstratified growth rate. The difference in <r max 
probably results from the constrasting global structure of 
the channel flows. While an unstratified parasite attacks 
an infinite lattice of channels of equal amplitude, its strat- 
ified cousins can only attack a finite number of channels, 
of varying amplitude. Consequently, the amount and config- 
uration of the free energy available are dissimilar, and ob- 
viously influence the maximum growth rates that parasites 
can achieve. 

Because the fastest growing modes have large k, the 
eigenfunctions tend to localise near the midplane where the 
velocity shear is greatest. They also lose any large-scale en- 
velope they may have exhibited. This can be observed in 
Fig. 12 which presents the fastest growing mode on the 
n = 21 channel when /3 — 1000. Here the mode is indi- 
vidually kinking the middle 12 jets, and ignoring the other 
8 at larger altitudes. Because these central jets are more 
closely packed they offer greater velocity shear for the para- 
sitic mode. This behaviour is typical for the fastest growing 
modes. It follows that, if they can be adequately resolved, 
parasites will be found primarily near the midplane in nu- 
merical simulations. 

A final point concerns the consistency of the interme- 
diate limit we employ. The (J max of the stratified parasites 
are not only smaller than their unstratified counterparts, 
they are sufficiently small to endanger the original scaling 
a ~ bQ. This means one should revise the lower limit in 
Eq. (23), which defines the intermediate regime, possibly 
increasing it to 10 or even 100. But unless /3 is very large 
indeed (> 10 , say), there may be no intermediate regime to 
speak of: one must either account for the background shear 
and channel growth, or one must include compressibility. 
Having said that, the general features of the analysis — the 
estimates of the growth rates and the modes' general struc- 
turem — should offer a satisfactory guide to the parasitic 
physics at any b. In the next subsection we discuss the more 
realistic non-intermediate regimes. 

4.5 Discussion 

Our treatment of the stratified parasitic modes provides 
only a limited survey of this problem's rich assortment of 
mode behaviour. But the salient dynamical points remain 
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Figure 12. Selected eigenfunction components of the kink mode 
upon the n = 21 channel for /3 = 1000. The horizontal wavcnum- 
ber is k = 0.575 K (= 12.36) and 9^ = 9, thus the growth rate is 
cr/(bQ) = 0.0329. The black vertical arrows indicate the locations 
of the MRI channels in the range 1.2 > z > 0, a symmetric set 
occurs for z < but is omitted (there are 20 altogether). Note 
that the mode has localised upon those channels near the mid- 
plane, while ignoring those at larger altitudes. This is typical of 
larger k modes, and hence the fastest growing parasites. 



clear: the fastest growing mode is always the hydrodynami- 
cal kink mode (with 9^ =9), favouring a wavenumber some 
half that of its host channel K. This is also the case in the 
unstratified analysis. While there exist an interesting array 
of other modes, exhibiting complicated morphologies, these 
pinch and kink-pinch modes grow significantly slower, and 
so should be dynamically subdominant. In short, if a chan- 
nel is to be disrupted by a parasite, it will be the kink mode 
doing the disrupting. Finally, the maximum growth rates 
that stratified parasitic modes attain are significantly less 
than those calculated from unstratified models. 



4-5.1 Small amplitude channels 

Our mathematical characterisation of the parasitic modes 
only holds in an intermediate limit of channel amplitude: 
channels must be sufficiently large that we may omit the 
Keplerian shear and the channels' exponential growth, but 
sufficiently small that compressible effects are negligible. 
However, in order to understand the role of the parasites 
in the MRI dynamics more generally, we need to extrapo- 
late these intermediate results into the regimes of small and 
large channel amplitudes. 

An important dynamical question is: How large will an 
MRI channel become before it is destroyed by a parasite? 
Will a channel reach equipartition before collapsing? This is 
important not only to understand (and possibly control) the 
MRI's initial eruptive behaviour, but also to assess whether 
continuous channel production and destruction are key in- 
gredients in MRI-induced turbulence (see Pessah and Good- 
man 2009, Pessah 2009). In order to answer this question we 



must understand the parasites' behaviour when the channel 
amplitudes are small. 

One could, of course, assume that the intermediate am- 
plitude profiles and growth rates (previously calculated) 
offer a reasonable approximation in the small amplitude 
regime, as do Pessah and Goodman (2009) and Pessah 
(2009). This means one expects the exponential growth of 
the channel and the Keplerian shear to not qualitatively al- 
ter the solutions. By simply equating the channel growth 
rate and the kink mode growth rate, we can obtain a lower 
bound on the minimum b necessary for channel disruption. 
In an unstratified ideal MHD channel this value is roughly 
4 (Pessah and Goodman 2009); in the stratified model em- 
ployed in the previous section with /3 = 1000, a similar cal- 
culation gives 6 m i n « 15. We now argue that both lower 
bounds are far too small. 

In LLB09 the influence of the Keplerian shear and chan- 
nel growth is discussed in a qualitative way. In that paper 
the channel growth was emphasised as key in delaying the 
onset of the parasitic modes. The role of the shear in doing 
the same was not considered. In fact, the shear should im- 
pede the parasites equally efficiently. The Keplerian shear 
introduces a time-dependence to the parasitic wavevector k. 
As time progresses the orientation angle 9k will decrease and 
the wavenumber k will increase. But parasitic modes only 
grow when their 8k lie within some range encompassing 6 
and when k < K (see Eq. (30) and Fig. 7, for example). It 
follows directly that these modes will grow only transiently, 
enjoying a short burst as 0k{t) and k(t) pass through their 
respective neighbourhoods of permitted growth. If the mode 
does not grow sufficiently fast within this finite time interval 
— sufficiently fast to reach an amplitude comparable to 6 — 
then the channel remains undisrupted. This time restriction 
provides a means to calculate the minimum 6 necessary for 
channel disruption. We present these details in Appendix 
B. There we find that for the unstratified problem, and for 
an initial parasitic amplitude 1/10 or 1/100 of the channel, 
the critical amplitude is b m i n « 24 or 40. In the stratified 
problem, these values are of order 100. Both estimates are 
considerably larger than if the Keplerian shear is omitted. 
Both also show that channels, particularly stratified chan- 
nels, are destined to achieve large amplitudes and will likely 
hit equipartition before the parasites can topple them: a 
channel that has grown from b — 1 to a hundred times that, 
will have decreased the midplane j3 by a factor 10 4 . 

4-5.2 Compressible channels 

Once the system has reached equipartition the analyses of 
Sections 2 and 4 break down. The profile of the MRI chan- 
nel will alter, its magnetic pressure sufficiently strong to 
squeeze mass into its jet centres (see Figs 3 and 5 and Stone 
and Miller 2000). The channel jets will become thinner and 
more intense, while simultaneously coinciding with narrow 
current sheets. If sufficiently thin, the influence of resistivity 
will become important, giving rise to new parasitic tearing 
modes (explored in LLB09). Finite resolution will ensure 
that the magnetic reconnection induced by these modes fea- 
tures prominently in stratified box simulations. Compress- 
ibility may introduce parasitic Parker modes, which will ex- 
ploit the strong bands of magnetic field between jet centres. 
The character of the original kink modes will change as well, 



though perhaps not drastically. As shown in LLB09, when 
the jets narrow and become more intense, the kink mode 
tends to grow on shorter horizontal length scales, and can 
split into multiple modes each localised to a single jet. This 
new wealth of dynamical behaviour should be observable in 
resolved 3D stratified simulations, and we intend to present 
some of this work in a later paper. 



5 CONCLUSION 

In summary, we have shown that channel flows, a subset 
of the linear axisymmetric MRI modes, are approximate 
nonlinear solutions in the vertically-stratified shearing box, 
when in the limit of large midplane /3. This generalises the 
better known unstratified result (GX94), but only holds 
when the equilibrium consists of a net vertical magnetic 
flux. We compute stratified channel flows numerically for 
an isothermal model, and confirm the basic result with ID 
simulations. The MRI eigenfunctions and growth rates, con- 
versely, offer a valuable numerical benchmark on the perfor- 
mance of codes operating in the vertically-stratified shearing 
box. In addition, these simulations provide estimates (for fu- 
ture 3D simulations) of the resolution required to resolve the 
fastest growing channels when using second order numerical 
schemes: 25 cells per scaleheight are needed when /3 = 100, 
50 cells when /3 = 10 3 , and more than 200 cells per scale 
height are mandatory when /3 = 10 4 . 

Stratified channel flows, like their unstratified counter- 
parts, are subject to a variety of magnetised shear flow in- 
stabilities, or 'parasitic modes'. Some of these we directly 
compute in an intermediate limit of channel amplitude. As 
in the unstratified setting, the fastest growing parasite is the 
'unmagnetised' kink mode. But it is unlikely to destroy its 
host before equipartition fields have been reached. 

At present we are extending these findings with 3D nu- 
merical simulations in stratified boxes, where the parasitic 
modes' behaviour can be confirmed, not only in the initial 
stages of a simulation but once it has relaxed into a satu- 
rated turbulent state. Because of the inevitably of large mag- 
netic fields, and also finite resolution, these simulations will 
permit a set of interesting physical effects which we have ne- 
glected here but which deserve further study. These include, 
most importantly, magnetic buoyancy and magnetic recon- 
nection. Also, these simulations will allow us to test the re- 
current formation of channels once the MRI has saturated in 
turbulence (Suzuki and Inutsuka 2009). In particular, they 
can establish whether larger radial domains impede recur- 
rent channel formation, by analogy with unstratified boxes 
(Bodo et al. 2008, LLB09). 

Lastly, we speculate on the potential role dominant 
channel flows play in real magnetised accretion disks. Global 
disk models support analogues of the stratified MRI chan- 
nels, localised in radius to narrow rings. Such modes might 
also exhibit a weaker version of the nonlinear property ex- 
amined in this paper, permitting one mode to dominate the 
initial evolution of small amplitude initial conditions. In- 
deed, simulations with a poloidal field piercing the disk re- 
port such behaviour (Hawley 2000, 2001). In fully-developed 
turbulent disks, on the other hand, the dominance of a single 
channel will be rare. But occasional 'run-away' channel flows 
might occur in the marginally stable fluid that divides re- 



gions of stable and unstable MRI — dead-zone interfaces in 
protostellar disks, for instance. It may take only a small fluc- 
tuation in the ionisation and thermal properties of the gas to 
ignite instability in this previously quiescent critical layer, 
and hence launch a powerful MRI channel flow. Disruptive 
events of this kind may control the integral-scale intermit- 
tency that should characterise the dead-zone interface, and 
may even be related to the initial phases of outburst be- 
haviour (Zhu et al. 2009). Similarly, channel flows may also 
feature in the onset of episodic accretion in dwarf nova disks 
(Gammie and Menou 1998). 
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APPENDIX A: ASYMPTOTIC ANALYSIS OF 
PARASITIC INSTABILITIES 

In this appendix we sketch out an asymptotic solution to Eqs 
(25)-(27) in the limit of small horizontal parasitic wavenum- 
ber: < k <C 1. We concentrate on a description of the kink 
mode as it is the fastest growing, and hence most dynami- 
cally significant. 

First we manipulate the governing equations into single 
second order equation for the pressure pertubation 

g{z) d(ld^_\ k ^ = (A1) 
yy ' dz \g(z) dz J y ' 

where 

g(z) = h(z) [a + ikMF n (z) cos A] 2 + |fc 2 G n (z) 2 sin 2 A. 

(A2) 

recall A = 6 — 6 k , which is just the difference in orienta- 
tion angles of the MRI channel and its parasite, and M is 
the Mach number. The algebra required to derive this can 
be greatly reduced by using Lagrangian displacements, not 
momenta (see GX94). The reader will notice that Eq (Al) is 
similar to Eq. (25) in GX94. We can reproduce this equation 
if we solve instead for the Lagrangian vertical displacement 
and assume the disk is unstratified. One can derive an anal- 
ogous 'semicircle theorem' for (Al) which shows that both 
the growth rate and frequency of the parasite is bounded 
above by k. Thus as k goes to 0, so must both the real and 
imaginary part of a. 

Al Kink modes 

We now suppose that k is sufficiently small that we may drop 
the last term in Eq. (Al). This assumption stipulates that 
the horizontal variation of the parasite is much greater than 



its vertical variation, as well as variations in the channel and 
the equilibrium density (which are of order H at most). 
We can solve directly, in this case, for 



* = ci + c 2 



g(z) dz 



where ci and C2 are integration constants and zb is the 
upper boundary of the disk (and may be oo). We require 
that = at 2 = ±zb (the vanishing boundary conditions). 
This furnishes ci = and the eigenvalue equation 



/* 



gdz — 0. 



After some manipulation we obtain the growth rate explic- 
itly: 



a 2 = A n k 2 M 2 [cos 2 A -MX 2 sin 2 A] , (A3) 



where 



= ( / ' /// ,' dz) I (J"" //-/;) . ( \ ! ) 



The Alfvenic Mach number Ma = uo/v A nas been intro- 
duced in the above, which can be re-expressed in terms of 
M and /3. In deriving equation (A3) the following two iden- 
tities were required 



f B hF 2 dz = f 

" — 2fi J — % 



G 2 dz, 



f 

J — ; 



hFdz = 0, 



the first follows from multiplying Eq. (17) by F and inte- 
grating, the second from the relationship between F and 
dG/dz. 

As anticipated, the growth rate scales like k and will be 
real (and positive) unless 



\9 - 6k\ > arctanMA, 



(A5) 



which is precisely the stability condition we find in the un- 
stratified case for kink modes (see GX94, and LLB09). The 
criterion just says that if a parasite is to grow, its wavevector 
must point sufficiently along the shear flow. 

The cigenfunctions can be obtained analytically if wc 
use the 'approximate isothermal model' (cf. Section 2.4.3). 
Then A n = l/(2n + 1), and in the hydrodynamical case of 
8k = 9 we obtain 

oc scch 2 z ^tanhz — iVtij , (A6) 

*2fc oc tanhzsech 4 z |^(5 — 4cosh2z) + 2iV5(l + cosh2z)j 

(A7) 

for n = 1 and n = 2 respectively. The eigenfunction is 
plotted in Fig. 8, and provides a check on the full numerical 
solution. 

Finally note that if the disk is unstratified we can com- 
pute A n = 1/2 if the limit of zb — > oo is taken carefully in 
both the denominator and numerator of (A4) (see also GX). 
Interestingly, this prefactor is independent of the mode num- 
ber or K, in contrast to the stratified case of the previous 
paragraph where larger n furnishes slower growing modes at 
small k. According to Section 4.5.1 this seems also to be the 
case for general k. 



A2 Pinch and kink-pinch modes 

The mathematical derivation of the pinch and kink-pinch 
modes is a great deal more involved and subtle, not least 
because the real part of the growth rate a comes in at higher 
orders. We will not produce the full analysis here. Instead, 
we work out the problem for the simple n = 2 case to leading 
order. 

The n — 2 MRI mode is a single jet centred at z = 0. 
It follows that there are two types of parasite: odd (kink 
modes) for which *(0) = 0, and even (pinch modes) for 
which d\l/ /dz(0) = 0. At leading order the latter boundary 
condition gives the equation: 

3(0) = 0, 

which is satisfied only if 

a = -ikMF(0) cos A. (A8) 

Applying the decaying boundary condition at z — zb com- 
pletes the description of the pinch mode, the vertical profile 
of which is 



* oc 1 



-a 



g{z)dz I 



f 

Jo 



gdz 



(A9) 



on the domain z £ [0, zb]. If we apply the approximate 
isothermal model and furthermore take 9 — 9k, then this 
reduces to 



^2r> oc 1 — tanh z. 



(A10) 



The eigenfunction is plotted in Fig. 9 against the numerical 
solution. The real part of the growth rate comes in at higher 
order and requires an additional asymptotic matching pro- 
cedure, the details of which we do not provide here. 



APPENDIX B: 
CHANNELS 



SMALL AMPLITUDE 



The purpose of this appendix is to derive an estimate for the 
smallest amplitude b that permits parasitic modes to over- 
come the primary channel flow. Channels with amplitudes 
below this critical value host parasites that grow too slowly 
to cause disruption. We define 'channel disruption' when 
the amplitudes of the channel and the parasite are equal, in 
contrast to when their growth rates are equal (Pessah and 
Goodman 2009, Pessah 2009). In the following, the influence 
of the background shear will be included, and we will find 
that it is a primary agent delaying channel disruption. We 
do not treat the influence of the Coriolis force. 

For simplicity, the disk is assumed to be unstratified. 
Consider the fastest growing MRI channel, with s — (3/4)0 
and 9 — n/4. Its fastest growing parasite, the one in which 
we will be most interested, is the kink mode. It achieves its 
maximum growth when 9k = 8 = tv/4 and k w 0.6 K, at 
which point a w 0.2 fcfl It, however, can grow at a reduced 
rate within the sector \9 — 9 k \ < # C rit < 7r/4, and for k < K 
(see GX94). 

Suppose that we allow for the background shear's ef- 
fect on the parasitic wavevector k. Under its influence the 
wavevector will become time-dependent, lengthening with 
time and rotating so that it points more and more radially. 
These changes are summarised by 



k x {t) — k x + — Qkyt, 



(Bl) 



where (k x , k y ) is the initial value of the wavevector. Note 10 4 . By then it is more than likely that equipartition has 
that the azimuthal component ky is independent of time. been reached. 

We analyse a fiducial kink mode with k y — (0.6/V2)K 
and k° x = 0. 1 Thus the kink mode starts out with a purely 
azimuthal wavevector and zero growth rate. It, however, will 
achieve the maximum growth possible at some later time 
when k x = k y , i.e. when k is perfectly aligned with the 
channel. After passing this 'sweet spot' its growth rate will 
decline, and when k x = \/0.82 % it will be zero, because at 
that point k — K. If we permit the mode to grow between 
orientations k x = and k x = k x , we can then estimate the 
total time the mode grows from Eq. (Bl). We have, in fact, 

The growth time alotted to a kink mode is thus roughly a 
quarter of an orbit. This is the maximum time available for 
it to disrupt its host. If it cannot reach a sufficiently large 
amplitude within this time then the channel will keep grow- 
ing unimpeded (until another kink mode emerges and tries 
its luck). Of course, the success of the parasite's attack de- 
pends on the average size of a over the period of growth. We 
can, in fact, easily derive a minimum value for this average 
growth rate, below which a kink mode cannot disrupt the 
channel. And because a is proportional to b, this in turn 
allows us to calculate the minimum 6 below which channels 
are safe from kink modes. 

To ease the calculation we shall take this average value 
to be half the maximum growth, i.e. cr av = 0.1 bfl. We also 
must attribute a starting amplitude to the channel, which 
we take to be Rb, where R is some fraction. During the time 
t = 1.42/n, the channel will grow by approximately e 1 ' 07 , 
while the kink mode will grow by approximately e°' 142i> . The 
ratio of their amplitudes at the end of this time is hence 

parasite amplitude ^ Q ^ & _ 

channel amplitude 

If this ratio is set to 1, then we can calculate the critical 
(minimum) b. We have simply 

6 cr « 7.54 - 7.04 In R. (B4) 

If the initial parasitic amplitude is set to one tenth the chan- 
nel, i.e. R = 1/10 then the minimum b for channel disruption 
is about 24, while R — 1/100 gives 6 cr closer to 40. In either 
case, channels will grow to significantly larger amplitudes 
than predicted by Pessah and Goodman (2009) and Pessah 
(2009), who neglect the influence of the background shear 
in retarding the parasitic growth. 

The analysis above was carried out for the unstratified 
problem, but it is easy to generalise to the stratified disk. 
As mentioned in Section 4, the stratified parasitic growth 
rates are smaller than unstratified ones. If this is taken into 
account we obtain b CI ~ 100 for both R = 1/10 and R = 
1/100. Consequently, by the time a parasite can make a 
successful attack, the channel will have reached enormous 
amplitudes, decreasing the midplane plasma beta by a factor 



1 We assume that a kink mode poses the greatest threat to chan- 
nel stability. Though pinch and kink-pinch modes grow for a 
larger range of 9k, they grow for a smaller range of k. This com- 
bined with their smaller growth rates leads us to discount them. 



